{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 133,
   "id": "1e2fcd89-41c4-43c3-9fce-2aaa13173349",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABHgAAAEYCAYAAAAnPkG+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAsOklEQVR4nO3dcZBc9XXg++/RSCJjQlY4CIwGHMhaO6+csEFrFw5FZXdsxwZrnUWhkl14r2KcuEp2XqgK+xLWaLcqTuxXa+1SyW52sSFam8KucgDXRsjaRYvgGffDrsIBg2QLGQZkwtqa4aGH7QEG5lnS6Lw/ukduDT0zrem+3ff2fD9VU9197+/ePn26+95fn7n3/iIzkSRJkiRJUnWt6ncAkiRJkiRJ6owFHkmSJEmSpIqzwCNJkiRJklRxFngkSZIkSZIqzgKPJEmSJElSxVngkSRJkiRJqjgLPJIkaUWLiAsj4msR8VREHIyIP2xMf3NEPBgRzzZuz15g+asiYjwiDkXEzb2NXpIkqS4ys98xSJIk9U1EnA+cn5lPRMRZwOPAFuDDwI8yc3ujcHN2Zn583rJDwDPA+4DDwGPAdZn53R6+BEmSJI/gkSRJK1tmvpCZTzTuvwo8BYwAVwNfaDT7AvWiz3yXAYcy87nMPArc3VhOkiSpp1b3O4BW1q1bl29729v6HcbAee211zjzzDP7HcZAMrfFMbfFMbfFMbfFePzxx1/KzPVFPkdEXARsAv4WOC8zX4B6ESgizm2xyAjwg6bHh4F3tVjvVmArwM/8zM+8461vfWuXIxfAiRMnWLXK/18WwdwWx9wWw7wWx9wW55lnnumor1PKAs95553Ht771rX6HMXBqtRpjY2P9DmMgmdvimNvimNvimNtiRMT/LHj9Pwv8DXBjZr4SEW0t1mLaG85/z8wdwA6A0dHRHB8f7yRULcDvXnHMbXHMbTHMa3HMbXE67etYdpMkSSteRKyhXtz5UmbubEx+sXF9nrnr9Bxpsehh4MKmxxcAk0XGKkmS1IoFHkmStKJF/VCdzwNPZeZfNM3aDVzfuH898JUWiz8GbIyIiyNiLXBtYzlJkqSessAjSZJWuiuA3wHeExH7G3+bge3A+yLiWeqjZG0HiIgNEbEHIDOPAzcAe6lfnPnLmXmwHy9CkiStbKW8Bo8kSVKvZOY3aH0tHYD3tmg/CWxuerwH2FNMdJIkSe3xCB5JkiRJkqSKs8AjSZIkSZJUcUsWeCLiwoj4WkQ8FREHI+IPG9PfHBEPRsSzjduzF1j+qogYj4hDEXFzt1+AJEmSJEnSStfONXiOA3+UmU9ExFnA4xHxIPBh4KuZub1RuLkZ+HjzghExBHyG+oUJDwOPRcTuzPxuN19EWe3aN8Ete8eZnJphw7phbrpylC2bRvodliRpwLn/kSRJWnmWPIInM1/IzCca91+lPkLECHA18IVGsy8AW1osfhlwKDOfy8yjwN2N5Qbern0TbNt5gImpGRKYmJph284D7No30e/QJEkDzP2PJEnSynRao2hFxEXAJuBvgfMy8wWoF4Ei4twWi4wAP2h6fBh41wLr3gpsBVi/fj21Wu10QiudT9VeZ+ZYnjJt5tgsn/rKt1n38rN9iWl6erryeS0rc1scc1scc1ucfua2jPsfSZIkFa/tAk9E/CzwN8CNmflKxEKjiZ66WItp2WIambkD2AEwOjqaY2Nj7YZWSj+6/77W0/+/pF+vrVar9e25B525LY65LY65LU4/c1vG/Y8kSZKK19YoWhGxhnpx50uZubMx+cWIOL8x/3zgSItFDwMXNj2+AJhcfrjVsWHd8GlNlySpG9z/SJIkrUztjKIVwOeBpzLzL5pm7Qaub9y/HvhKi8UfAzZGxMURsRa4trHcwLvpylGG1wydMm14zRA3XTnap4gkSSuB+x9JkqSVqZ0jeK4Afgd4T0Tsb/xtBrYD74uIZ6mPkrUdICI2RMQegMw8DtwA7KV+ceYvZ+bBAl5H6WzZNMKnr7mEtUP1FI+sG+bT11ziKCaSpEK5/5EkSVqZlrwGT2Z+g9bX0gF4b4v2k8Dmpsd7gD3LDbDKtmwa4a5Hvw/APR+9vM/RSJJWCvc/kiRJK09b1+CRJEmSJElSeVngkSRJkiRJqjgLPJIkSZIkSRVngUeSJEmSJKniLPBIkiRJkiRVnAUeSZIkSZKkirPAI0mSJEmSVHEWeCRJkiRJkirOAo8kSZIkSVLFWeCRJEmSJEmqOAs8kiRJkiRJFbe63wFIkiT1U0TcAXwQOJKZv9yYdg8w2miyDpjKzEtbLPs88CowCxzPzHf2IGRJkqQ3sMAjSZJWujuBW4Evzk3IzH8xdz8i/hx4eZHl352ZLxUWnSRJUhss8EiSpBUtMx+OiItazYuIAP458J6eBiVJknSaLPBIkiQt7NeAFzPz2QXmJ/BARCTwV5m5o1WjiNgKbAVYv349tVqtiFhXvOnpaXNbEHNbHHNbDPNaHHNbXhZ4JEmSFnYdcNci86/IzMmIOBd4MCKezsyH5zdqFH52AIyOjubY2Fghwa50tVoNc1sMc1scc1sM81occ1tejqIlSZLUQkSsBq4B7lmoTWZONm6PAPcCl/UmOkmSpFNZ4JEkSWrt14GnM/Nwq5kRcWZEnDV3H3g/8GQP45MkSTrJAo8kSVrRIuIu4BFgNCIOR8RHGrOuZd7pWRGxISL2NB6eB3wjIr4NPArcl5n39ypuSZKkZktegyci7gA+CBzJzF9uTLsHGG00WQdMZealLZZ9HngVmAWOZ+Y7uxK1JElSl2TmdQtM/3CLaZPA5sb954BfKTQ4SZKkNrVzkeU7gVuBL85NyMx/MXc/Iv4ceHmR5d+dmS8tN0BJkiRJkiQtbskCT2Y+HBEXtZoXEQH8c+A9XY5LkiRJkiRJbep0mPRfA17MzGcXmJ/AAxGRwF81hghtKSK2AlsB1q9fT61W6zC0cpiamgEoxeuZnp4uRRyDyNwWx9wWx9wWpwy5LdP+R5IkScXrtMBzHfMuPjjPFZk5GRHnAg9GxNOZ+XCrho3izw6A0dHRHBsb6zC0crht/BEAxsYu73Mk9U7+oOS1bMxtccxtccxtccqQ2zLtfyRJklS8ZY+iFRGrgWuAexZq07gQIZl5BLgXuGy5zydJkiRJkqTWOhkm/deBpzPzcKuZEXFmRJw1dx94P/BkB88nSZIkSZKkFpYs8ETEXcAjwGhEHI6IjzRmXcu807MiYkNE7Gk8PA/4RkR8G3gUuC8z7+9e6JIkSZIkSYL2RtG6boHpH24xbRLY3Lj/HPArHcYnSZIkSZKkJXRyipYkSZIkSZJKwAKPJEmSJElSxVngkSRJkiRJqjgLPJIkSZIkSRVngUeSJEmSJKniLPBIkiRJkiRVnAUeSZIkSZKkirPAI0mSJEmSVHEWeCRJkiRJkirOAo8kSZIkSVLFWeCRJEmSJEmqOAs8kiRJkiRJFWeBR5IkSZIkqeIs8EiSJEmSJFWcBR5JkrSiRcQdEXEkIp5smvanETEREfsbf5sXWPaqiBiPiEMRcXPvopYkSTqVBR5JkrTS3Qlc1WL6f8jMSxt/e+bPjIgh4DPAB4C3A9dFxNsLjVSSJGkBFngkSdKKlpkPAz9axqKXAYcy87nMPArcDVzd1eAkSZLatLrfAUiSJJXUDRHxIeBbwB9l5o/nzR8BftD0+DDwrlYrioitwFaA9evXU6vVuh+tmJ6eNrcFMbfFMbfFMK/FMbflZYFHkiTpjW4DPgVk4/bPgd+b1yZaLJetVpaZO4AdAKOjozk2Nta1QPVTtVoNc1sMc1scc1sM81occ1teS56i5YUHJUnSSpOZL2bmbGaeAP4L9dOx5jsMXNj0+AJgshfxSZIkzdfONXjuxAsPSpKkFSQizm96+JvAky2aPQZsjIiLI2ItcC2wuxfxSZIkzbdkgccLD0qSpEEWEXcBjwCjEXE4Ij4C/PuIOBAR3wHeDfzLRtsNEbEHIDOPAzcAe4GngC9n5sG+vAhJkrTidXINnq5deBAG9+KDU1MzAKV4PV4MqzjmtjjmtjjmtjhlyG2Z9j9ll5nXtZj8+QXaTgKbmx7vAd5wJLMkSVKvLbfA09ULD8LgXnzwtvFHABgbu7zPkXgxrCKZ2+KY2+KY2+KUIbdl2v9IkiSpeO1cg+cNvPCgJEmSJElSeSyrwOOFByVJkiRJkspjyVO0GhceHAPOiYjDwCeAsYi4lPopV88DH2203QB8LjM3Z+bxiJi78OAQcIcXHpQkSZIkSeq+JQs8XnhQkiRJkiSp3JZ1ipYkSZIkSZLKwwKPJEmSJElSxVngkSRJkiRJqjgLPJIkSZIkSRVngUeSJEmSJKniLPBIkiRJkiRVnAUeSZIkSZKkirPAI0mSJEmSVHEWeCRJkiRJkirOAo8kSZIkSVLFWeCRJEmSJEmqOAs8kiRJkiRJFWeBR5IkSZIkqeIs8EiSJEmSJFWcBR5JkiRJkqSKs8AjSZIkSZJUcRZ4JEnSihYRd0TEkYh4smnaLRHxdER8JyLujYh1Cyz7fEQciIj9EfGtngUtSZI0jwUeSZK00t0JXDVv2oPAL2fmPwSeAbYtsvy7M/PSzHxnQfFJkiQtyQKPJEla0TLzYeBH86Y9kJnHGw+/CVzQ88AkSZJOw+qlGkTEHcAHgSOZ+cuNabcAvwEcBb4H/G5mTrVY9nngVWAWOO5/tiRJUgX9HnDPAvMSeCAiEvirzNzRqlFEbAW2Aqxfv55arVZEnCve9PS0uS2IuS2OuS2GeS2OuS2vJQs81A9bvhX4YtO0B4FtmXk8Iv4d9cOWP77A8u/OzJc6ilKSJKkPIuLfAMeBLy3Q5IrMnIyIc4EHI+LpxhFBp2gUfnYAjI6O5tjYWFEhr2i1Wg1zWwxzWxxzWwzzWhxzW15LnqLlYcuSJGkliojrqR/F/L9lZrZqk5mTjdsjwL3AZb2LUJIk6afaOYJnKR0ftgyDe+jy1NQMQClej4fSFcfcFsfcFsfcFqcMuS3T/qeKIuIq6kcn/5PMfH2BNmcCqzLz1cb99wOf7GGYkiRJJ3VU4OnWYcswuIcu3zb+CABjY5f3ORIPpSuSuS2OuS2OuS1OGXJbpv1P2UXEXcAYcE5EHAY+Qf308zOo918AvpmZH4uIDcDnMnMzcB5wb2P+auCvM/P+PrwESZKk5Rd4mg5bfm87hy1HxNxhyy0LPJIkSf2Qmde1mPz5BdpOApsb958DfqXA0CRJktq2rGHSmw5b/meLHbYcEWfN3ad+2PKTyw1UkiRJkiRJrS1Z4GkctvwIMBoRhyPiI9RH1TqL+mHL+yPi9kbbDRGxp7HoecA3IuLbwKPAfR62LEmSJEmS1H1LnqLlYcuSJEmSJEnltqxTtCRJkiRJklQeFngkSZIkSZIqzgKPJEmSJElSxVngkSRJkiRJqjgLPJIkSZIkSRVngUeSJEmSJKniLPBIkiRJkiRVnAUeSZIkSZKkirPAI0mSJEmSVHEWeCRJkiRJkirOAo8kSZIkSVLFWeCRJEmSJEmqOAs8kiRJkiRJFWeBR5IkSZIkqeIs8EiSJEmSJFWcBR5JkiRJkqSKs8AjSZIkSZJUcRZ4JEmSJEmSKs4CjyRJkiRJUsUtWeCJiDsi4khEPNk07c0R8WBEPNu4PXuBZa+KiPGIOBQRN3czcLVv174Jrtj+EB++/zWu2P4Qu/ZN9DskSRpIbm+ryb5O9fndk6Tiua0tv3aO4LkTuGretJuBr2bmRuCrjceniIgh4DPAB4C3A9dFxNs7ilanbde+CbbtPMDE1AwAE1MzbNt5wC+jJHWZ29tKuxP7OpXld0+Siue2thqWLPBk5sPAj+ZNvhr4QuP+F4AtLRa9DDiUmc9l5lHg7sZy6qFb9o4zc2z2lGkzx2a5Ze94nyKSpMHk9ra67OtUm989SSqe29pqWL3M5c7LzBcAMvOFiDi3RZsR4AdNjw8D71pohRGxFdgKsH79emq12jJDK5epRoWzX69nrsLaavqg5LgMpqenzWdBzG1xzG13lW172+/9zwDoal9nUPs5ZVC2796gcp9RHHNbDPPaXW5rq2G5BZ52RItpuVDjzNwB7AAYHR3NsbGxgsLqrdvGHwFgbOzyvjz/yDcfavllHFk3zKDkuAxqtZr5LIi5LY657a6ybW/7vf9ZIdru6wxqP6cMyvbdG1TuM4pjbothXrvLbW01LHcUrRcj4nyAxu2RFm0OAxc2Pb4AmFzm82mZbrpylOE1Q6dMG14zxE1XjvYpIkkaTG5vB459nYrwuydJxXNbWw3LLfDsBq5v3L8e+EqLNo8BGyPi4ohYC1zbWE49tGXTCJ++5hLWDtXf6pF1w3z6mkvYsmmkz5FJ0mBxeztw7OtUhN89SSqe29pqWPIUrYi4CxgDzomIw8AngO3AlyPiI8D3gd9utN0AfC4zN2fm8Yi4AdgLDAF3ZObBYl6GFrNl0wh3Pfp9pqam2Pvx9/Q7HEkaWG5vq8m+TvX53ZOk4rmtLb8lCzyZed0Cs97bou0ksLnp8R5gz7KjkyRJKph9HUmSNAiWe4qWJEmSJEmSSsICjyRJkiRJUsVZ4JEkSZIkSao4CzySJEmSJEkVZ4FHkiRJkiSp4izwSJIkSZIkVZwFHkmSJEmSpIqzwCNJkiRJklRxFngkSZIkSZIqzgKPJEmSJElSxVngkSRJkiRJqjgLPJIkSZIkSRVngUeSJEmSJKniLPBIkiRJkiRVnAUeSZIkSZKkirPAI0mSJEmSVHEWeCRJkiRJkirOAo8kSZIkSVLFWeCRJEmSJEmquGUXeCJiNCL2N/29EhE3zmszFhEvN7X5k44jliRJkiRJ0ilWL3fBzBwHLgWIiCFgAri3RdOvZ+YHl/s8kiRJkiRJWly3TtF6L/C9zPyfXVqfJEmSJEmS2rTsI3jmuRa4a4F5l0fEt4FJ4I8z82CrRhGxFdgKsH79emq1WpdC66+pqRmAvr+eqakZZmdn+x7HoJqenja3BTG3xTG3xSjL9rYs+58qi4hR4J6mSb8I/Elm/semNmPAV4C/a0zamZmf7FGIkiRJJ3Vc4ImItcA/A7a1mP0E8AuZOR0Rm4FdwMZW68nMHcAOgNHR0RwbG+s0tFK4bfwRAMbGLu97HFNTUwxKXsumVquZ24KY2+KY22KUZXtblv1PlXk6uiRJqpJunKL1AeCJzHxx/ozMfCUzpxv39wBrIuKcLjynJElSL3k6uiRJKrVunKJ1HQucnhURbwFezMyMiMuoF5R+2IXnlCRJ6qWOTkcf1FPRy6Qsp0cOKk/rLY65LYZ5LYbb2nLrqMATEW8C3gd8tGnaxwAy83bgt4Dfj4jjwAxwbWZmJ88pSZLUS904HX1QT0Uvk7KcHjmoPK23OOa2GOa1GG5ry62jAk9mvg78/LxptzfdvxW4tZPnkCRJ6rNFT0dvur8nIj4bEedk5ks9jVCSJK143RomXZIkaVAtejp6RETjvqejS5KkvunWMOmSJEkDx9PRJUlSVVjgkSRJWoCno0uSpKrwFC1JkiRJkqSKs8AjSZIkSZJUcRZ4JEmSJEmSKs4CjyRJkiRJUsVZ4JEkSZIkSao4CzySJEmSJEkVZ4FHkiRJkiSp4izwSJIkSZIkVZwFHkmSJEmSpIqzwCNJkiRJklRxFngkSZIkSZIqzgKPJEmSJElSxVngkSRJkiRJqjgLPJIkSZIkSRVngUeSJEmSJKniLPBIkiRJkiRV3OpOFo6I54FXgVngeGa+c978AP4S2Ay8Dnw4M5/o5DkltbZr3wS37B1ncmqGDeuGuenKUbZsGul3WGrw/ZEkSZJUpI4KPA3vzsyXFpj3AWBj4+9dwG2NW0ldtGvfBNt2HmDm2CwAE1MzbNt5AMAiQgn4/kiSJEkqWtGnaF0NfDHrvgmsi4jzC35OacW5Ze/4yeLBnJljs9yyd7xPEamZ748kSZKkonV6BE8CD0REAn+VmTvmzR8BftD0+HBj2gvzVxQRW4GtAOvXr6dWq3UYWjlMTc0A9P31TE3NMDs72/c4BtX09HRfczvR+Jy1ml7197zfue2Gsr4/g5DbMirL9rYs+x9JkiT1RqcFnisyczIizgUejIinM/PhpvnRYplstaJGcWgHwOjoaI6NjXUYWjncNv4IAGNjl/c9jqmpKQYlr2VTq9X6mtuRbz7Usogwsm648u95v3PbDWV9fwYht2VUlu1tWfY/kiRJ6o2OTtHKzMnG7RHgXuCyeU0OAxc2Pb4AmOzkOSW90U1XjjK8ZuiUacNrhrjpytE+RaRmvj9SdUXE8xFxICL2R8S3WsyPiPhPEXEoIr4TEf+oH3FKK8GufRNcsf0hLr75Pq7Y/hC79k30OyQ18f2R+m/ZR/BExJnAqsx8tXH//cAn5zXbDdwQEXdTv7jyy5n5htOzJHVm7kK9/+q/foejsycYcZSmUvH9kSrPASWkPnPAgnLz/ZHKoZNTtM4D7q2PhM5q4K8z8/6I+BhAZt4O7KE+RPoh6sOk/25n4UpayJZNI9z16PcBuOejnpJRNr4/0sA6OaAE8M2IWBcR5/sPLam7FhuwwAJC//n+SOWw7AJPZj4H/EqL6bc33U/gD5b7HJIkSX3WlQElBnUwiTIpywXOB1W/L8xf1gELuqHfue2GMr4/g5DXMnJbW26dXmRZkiRpkHVlQIlBHUyiTMpygfNB1e8L85d1wIJu6Hduu6GM788g5LWM3NaWW0cXWZYkSRpkDighlYMDFpSb749UDhZ4JEmSWoiIMyPirLn71AeUeHJes93Ahxqjaf0qDighFWLLphE+fc0lrB2q/3wZWTfMp6+5xOu7lITvj1QOnqIlSZLUmgNKSCXigAXl5vsj9Z8FHkmSpBYcUEKSJFWJp2hJkiRJkiRVnAUeSZIkSZKkirPAI0mSJEmSVHEWeCRJkiRJkirOAo8kSZIkSVLFWeCRJEmSJEmqOAs8kiRJkiRJFWeBR5IkSZIkqeJW9zsASSvPrn0T3LJ3nMmpGTasG+amK0fZsmmkb+uRJEmSpKqzwCOpp3btm2DbzgPMHJsFYGJqhm07DwCcVnGmW+uRJEmSpEHgKVqSeuqWveMnizJzZo7Ncsve8b6sR5IkSZIGgQUeST01OTVzWtOLXo8kSZIkDQILPJJ6asO64dOaXvR6JEmSJGkQLLvAExEXRsTXIuKpiDgYEX/Yos1YRLwcEfsbf3/SWbiSqu6mK0cZXjN0yrThNUPcdOVoX9YjSZIkSYOgk4ssHwf+KDOfiIizgMcj4sHM/O68dl/PzA928DySBsjcBZD/1X/9DkdnTzCyzNGvurUeSZIkSRoEyy7wZOYLwAuN+69GxFPACDC/wCNJp9iyaYS7Hv0+APd89PK+r0eSJEmSqq4rw6RHxEXAJuBvW8y+PCK+DUwCf5yZB7vxnJIkSZJURrv2TXDL3nEmp2bYsMyjjLuxDkkrS8cFnoj4WeBvgBsz85V5s58AfiEzpyNiM7AL2LjAerYCWwHWr19PrVbrNLRSmGqM6NPv1zM1NcPs7Gzf4xhU09PTpchtWT5v7Wg31qVyO4ivuVfK8rkdNGXZ3pbt8yZJK8WufRNs23mAmWOzAExMzbBt5wGAtgs03ViHpJWnowJPRKyhXtz5UmbunD+/ueCTmXsi4rMRcU5mvtSi7Q5gB8Do6GiOjY11Elpp3Db+CABjY/09feS28UeYmppiUPJaNrVarRS5LcvnrR3txrpUbgfxNfdKWT63g6Ys29uyfd4kaaW4Ze/4ycLMnJljs9yyd7zt4kw31iFp5elkFK0APg88lZl/sUCbtzTaERGXNZ7vh8t9TkmSJEkqs8nGEZTtTi9qHZJWnk6O4LkC+B3gQETsb0z718BbATLzduC3gN+PiOPADHBtZmYHzylJkiRJpbVh3TATLQoxG9YN93QdklaeTkbR+gYQS7S5Fbh1uc8hSZIkSVVy05Wjp1w/B2B4zRA3XTna03VIWnmWfYqWJEnSIIuICyPiaxHxVEQcjIg/bNFmLCJejoj9jb8/6Ueskspjy6YRPn3NJawdqv/UGlk3zKevueS0rp3TjXVIWnm6Mky6JEnSADoO/FFmPhERZwGPR8SDmfndee2+npkf7EN8kkpqy6YR7nr0+wDc89HlXey+G+uQtLJ4BI8kSVILmflCZj7RuP8q8BTgv88lSVIpeQSPJEnSEiLiImAT8LctZl8eEd8GJoE/zsyDLZbfCmwFWL9+PbVarbhgV6ipqRlmZ2fNbUGmp6dLkdupxoWHyxDLUtqNdbHcVun1QrniLctndtC4rS03CzySJEmLiIifBf4GuDEzX5k3+wngFzJzOiI2A7uAjfPXkZk7gB0Ao6OjOTY2VmjMK9Ft448wNTWFuS1GrVYrRW5vG38EgLGx8p+y1G6si+W2Sq8XyhVvWT6zg8Ztbbl5ipYkSdICImIN9eLOlzJz5/z5mflKZk437u8B1kTEOT0OU5IkyQKPJElSKxERwOeBpzLzLxZo85ZGOyLiMup9qx/2LkpJkqQ6T9GSJElq7Qrgd4ADEbG/Me1fA28FyMzbgd8Cfj8ijgMzwLWZmX2IVZIkrXAWeCRJklrIzG8AsUSbW4FbexORJEnSwjxFS5IkSZIkqeI8gqcAu/ZNcMvecSanZlgztIoLzx7ud0iSpBXopVd/whXbH2JyaoYN64a56cpRtmwa6XdYkiRJKoAFni7btW+CbTsPMHNsFoCjsyf43kuvcdHN9zFi51qSVLC5fzJMTM2cMn1iaoZtOw8AuB+SJEkaQJ6i1WW37B0/WdyZb65zvWvfRI+jkiStBHP/ZJhf3Jkzc2yWW/aO9zgqSZIk9YIFni6bXKBTPcfOtSSpKIv9k2HOUvspSZIkVZMFni7bsG7p6+3YuZYkFaGd/Us7+ylJkiRVjwWeLtq1b4LXfnJ8yXarIjxNS5LUVbv2TbAqFh3RG4DXjx53HyRJkjSAvMhyl8y/uPJiZjN7dqHL5ottBrDpkw8w9foxR1ORpC5qHj1xVUACF998X8+2tXP7oNnMJdv++PVjXmxZkiRpAFng6dCufRP86e6DTM0cazl/KGC2RX975tgsN96znxvv2V/Y6Frzi05JvWMPjqYiSd0yf1vbvM0vclu70GhZ861eFRw/ceqOaG4f9Gf/7SCf+I1fcj8gSZI0ACzwnIalijmttCruzDcxNXOy2DMUwWzmydvTLf602+GH3hSZJGkQ9Wpb2/w88/cP7Zpf3Gn249ePceM9+/k/vryfE8my9z2SJEnqv44KPBFxFfCXwBDwuczcPm9+NOZvBl4HPpyZTyy13udfOcFFN9/HquCUDmdQPwoFODmv1bTF2i93XvO007F2aBXHZk+0vexcp33utrn4027My9GqyNTPfFcyhr17+h7D2qFVrBtezaV/9sDJQmS/8t3OD8SXXv0JV2x/6JQfry1juP++BZ9n7dAqLjx78YvGtvqR3I/P3epVwZvftObkay7FZ79Fbkv7HStRDMvVvK1tJ4Zm8/cP7Vg7VL/U3tHZE4u2m3uuTvY9vXjP177lbe9o+8VLkiStMMsu8ETEEPAZ4H3AYeCxiNidmd9tavYBYGPj713AbY3btszvcGaLea2mLdZ+ufOW051fFXDh2cMcnnqdY7P5hs766Wo35k6VId/GsLwYjs6e4Mj0UWgxr9evdbEfsq1+JC83hqOzJ/jeS69x0c33vWH9rQqf/XzPj5/IU96fsn7uyvr5LlMM3dBODJ2aK37+3Q9fW/Y6y5Dv7mVdkiRpcHUyitZlwKHMfC4zjwJ3A1fPa3M18MWs+yawLiLO7+A5K+Xinz+Tc846g7Wrh7j458/sdzhS38z/sdbNH8nN5q/fH4Vayc5YHZxz1hmcc9YZ7oMkSZJWgE5O0RoBftD0+DBvPDqnVZsR4IXFVnzBq0f4t1//bAeh9d/wmiHWrK7Xz177yXHOPGM1x46fYOb4rL86JUnFCRhePcTR2ROcecZPd/PHjp9oa6THMvtQvwOQJEkqsU4KPNFi2vzSRTtt6g0jtgJbAX7xzJ/rIKz+GloF575pFX9vLUDjmgfDq07ef/noKl56PTnWrePv2xQB5w1DrFrF//PaCQo6gEKSVrQIeMuZq8gTJ3hxhp5va9esCs55UzT2QT/d98x5+egqjrx2oq0BACRJklQtnRR4DgMXNj2+AJhcRhsAMnMHsAPgjPM35sd/7X/vILTemLuuyHJHG+nG6CinE1++/Cz/ZGzstEZ/kSS1Z25b+483jVCr1Vj79zae3NZ2eiH8Zp2Otjin1T6o9J67sd8RSJIklVYnBZ7HgI0RcTEwAVwL/K/z2uwGboiIu6mfvvVyZi56elazMo600s2hY7dsGllwPZ2M+nP2m9bwid/4pTesu1Z79g3Pu1iRqQz5NobTi2Fk3TDv/l/W89+//ULfRtE63R+yc+27EUM7uvG6lpvTs9+0hn/6D8/na0//v+UZRauPz1PlGNrZJ8zfxu/aN8Gf7j64rO9mEcOWt9oHdWPEuaLfc0mSJLW27AJPZh6PiBuAvdSHSb8jMw9GxMca828H9lAfIv0Q9WHSf7eddV/0c6sY3/5PlxvaQFis+FPF5xl0tVqNsbGxfodx0v+55ZK+Pv9SP2RP50frUrltt0i5UOFzJSvb53bQVWF7W/YY49998PF+xyBJklRWnRzBQ2buoV7EaZ52e9P9BP6gk+eQVD29/JFY9h+kkiRJktQLnQyTLkmSJEmSpBKwwCNJkiRJklRxFngkSZIkSZIqrqNr8EiSJA2yiLgK+EvqA0p8LjO3z5sfjfmbqQ8o8eHMfGKxdT7/ygn+/rY9K270uaLnrV4VvGkoufTPHujaSJJlz0PP833/fX2PYSggIrjo5vt6/r6e7nuxdmgV64ZXt/eZvP++BQelGArY9MkH+PHrx0r7Whd7f/r62d973+B8//oYQ6tt7RXbH1rWqLBlf639jmHtW972DjpggUeSJKmFiBgCPgO8DzgMPBYRuzPzu03NPgBsbPy9C7itcbuo2ax357Jp2onsbF4311XFGI6fSF45ARw71rcYep2HlfiezyZw8n5vX+vp5uHo7AmOTB+FFvPae60/fc0/fv1YT1/Pct/zVu/PIHzuVnoMzdPmtrWvTM30LYZBz3cnLPBIkiS1dhlwKDOfA4iIu4GrgeYCz9XAFxsjh34zItZFxPmZ+cJCK73g1SP8269/tsi4JUlSBX2ow+VLWeB55plnpiNivN9xDKBzgJf6HcSAMrfFMbfFMbfFMbfFGO3x840AP2h6fJg3Hp3Tqs0IcEqBJyK2AlsBVg3/HB967pmuBytJkqrt+MtHOlq+lAUeYDwz39nvIAZNRHzLvBbD3BbH3BbH3BbH3BYjIr7V66dsMW3+AdTttCEzdwA7oP46fvL6y34+CuB3rzjmtjjmthjmtTjmtjid9nUcRUuSJKm1w8CFTY8vACaX0UaSJKlwFngkSZJaewzYGBEXR8Ra4Fpg97w2u4EPRd2vAi8vdv0dSZKkopT1FK0d/Q5gQJnX4pjb4pjb4pjb4pjbYvQ0r5l5PCJuAPZSHyb9jsw8GBEfa8y/HdhDfYj0Q9SHSf/dNlbt56M45rY45rY45rYY5rU45rY4HeU2MrswFpckSZIkSZL6xlO0JEmSJEmSKs4CjyRJkiRJUsX1tcATEb8dEQcj4kREvLNp+kURMRMR+xt/tzfNe0dEHIiIQxHxnyKi1fCkK95CuW3M29bI33hEXNk03dyepoj404iYaPqsbm6a1zLPak9EXNXI3aGIuLnf8VRdRDzf+H7vnxt+MSLeHBEPRsSzjduz+x1nFUTEHRFxJCKebJq2YC7dFrRvgdxWejtrX6cY9nN6o+rfv7Kzr9Nd9nW6x75OcYru6/T7CJ4ngWuAh1vM+15mXtr4+1jT9NuArcDGxt9VxYdZSS1zGxFvpz4KyC9Rz91nI2KoMdvcLs9/aPqs7oEl86wlNHL1GeADwNuB6xo5VWfe3ficzv0Yuhn4amZuBL7aeKyl3ckbt48tc+m24LTdSet9T5W3s/Z1imE/p3eq/P0rLfs6hbGv0x13Yl+nKHdSYF+nrwWezHwqM8fbbR8R5wM/l5mPZP3q0F8EthQVX5Utkturgbsz8yeZ+XfUR/24zNx2Xcs89zmmKrkMOJSZz2XmUeBu6jlVd10NfKFx/wv4nW9LZj4M/Gje5IVy6bbgNCyQ24VUIrf2dYphP6fvKvH9Kzn7Or1hX2cZ7OsUp+i+Tr+P4FnMxRGxLyL+74j4tca0EeBwU5vDjWlq3wjwg6bHczk0t8t3Q0R8p3G43dyhigvlWe0xf92XwAMR8XhEbG1MOy8zXwBo3J7bt+iqb6Fc+lnujkHdztrX6T77Od03qN+/fjOH3Wdfp1j2dYrVlW3t6qKimxMR/xfwlhaz/k1mfmWBxV4A3pqZP4yIdwC7IuKXgFbnSq/Ycd6XmduFcmhuF7BYnqkf7v0p6rn6FPDnwO9hPjtl/rrvisycjIhzgQcj4ul+B7RC+FnuXOm3s/Z1imE/pzfs5/SNOew++zr94We5c13b1hZe4MnMX1/GMj8BftK4/3hEfA/4B9QrVhc0Nb0AmOxGnFW0nNxSz+GFTY/ncmhuF9BuniPivwD/vfFwoTyrPeavyzJzsnF7JCLupX5454sRcX5mvtA4feFIX4OstoVy6We5Q5n54tz9sm5n7esUw35Ob9jP6Rtz2GX2dQpnX6cg3ezrlPIUrYhYP3fxoIj4ReoXwnuucSjYqxHxqxERwIeAhf6Do9Z2A9dGxBkRcTH13D5qbpensXGb85vUL/oIC+S51/FV2GPAxoi4OCLWUr+42O4+x1RZEXFmRJw1dx94P/XP6m7g+kaz6/E734mFcum2oEODup21r1MY+zldNKjfv5Kwr9NF9nV6wr5OQbq5rS38CJ7FRMRvAv8ZWA/cFxH7M/NK4B8Dn4yI48As8LHMnLsQ0e9Tv/L0MPA/Gn+aZ6HcZubBiPgy8F3gOPAHmTnbWMzcnr5/HxGXUj9U7nngowBL5FlLyMzjEXEDsBcYAu7IzIN9DqvKzgPurf+mYTXw15l5f0Q8Bnw5Ij4CfB/47T7GWBkRcRcwBpwTEYeBTwDbaZFLtwWnZ4HcjlV5O2tfpxj2c3rGfk5B7Ot0nX2dLrKvU5yi+zpRH0hAkiRJkiRJVVXKU7QkSZIkSZLUPgs8kiRJkiRJFWeBR5IkSZIkqeIs8EiSJEmSJFWcBR5JkiRJkqSKs8AjSZIkSZJUcRZ4JEmSJEmSKu7/B14w3XAVTUUgAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1152x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import math\n",
    "\n",
    "def getSin(amp,freq,phase,sampleList):\n",
    "    return amp*np.sin(-2*math.pi*freq*sampleList+phase)\n",
    "def getCos(amp,freq,phase,sampleList):\n",
    "    return amp*np.cos(-2*math.pi*freq*sampleList+phase)\n",
    "# 1. 获得混合波形\n",
    "srate=3000\n",
    "t=np.linspace(0,1,srate)\n",
    "s1=getSin(amp=5, freq=30, phase=0,sampleList=t)    \n",
    "s2=getCos(amp=3, freq=5, phase=0,sampleList=t)  \n",
    "s3=getSin(amp=10, freq=100, phase=0,sampleList=t)\n",
    "s4=getCos(amp=20, freq=120, phase=0,sampleList=t)  \n",
    "m=s1+s2+s3+s4\n",
    "# 2. 获得傅里叶系数\n",
    "fCoefs = np.fft.fft(m,srate)\n",
    "# 3. 获得振幅列表：每一个绕线的重心到原点的距离\n",
    "amp_list=2*np.abs(fCoefs/srate)\n",
    "# 把频率轴从0~1000 转变成 0~499 然后 -500~-1\n",
    "freqs = np.fft.fftfreq(len(amp_list), 1/srate)\n",
    "# 然后把 频率轴 和 数据 都变成 0hz 在中间，向左是负频率，向右是正频率的形式\n",
    "amp_shifted=np.fft.fftshift(amp_list)\n",
    "freq_shift=np.fft.fftshift(freqs)\n",
    "amp_shifted_copy = np.copy(amp_shifted)#获取图像\n",
    "amp_shifted_copy[amp_shifted < 0.5] = 0#使用掩模\n",
    "amp_shifted_copy[(freq_shift >= 110) | (freq_shift <= -110)] = 0\n",
    "# 把振幅列表画出来\n",
    "fig,ax=plt.subplots(1,2,figsize=(16,4))\n",
    "ax[0].grid()\n",
    "ax[0].stem(freq_shift,amp_shifted)\n",
    "ax[0].set_xlim([-150,150])\n",
    "ax[1].stem(freq_shift,amp_shifted_copy)\n",
    "ax[1].set_xlim([-150,150])\n",
    "ax[1].set_ylim([0, 20])\n",
    "ax[1].grid()\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "6f1506ce-30ef-4e2f-bb9d-187e54f99e22",
   "metadata": {},
   "outputs": [
    {
     "ename": "AttributeError",
     "evalue": "'numpy.ndarray' object has no attribute 'plot'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mAttributeError\u001b[0m                            Traceback (most recent call last)",
      "\u001b[0;32m/tmp/ipykernel_47/543279754.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m     34\u001b[0m \u001b[0mfig\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0max\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msubplots\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mfigsize\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m16\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m4\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     35\u001b[0m \u001b[0;31m#左上角子图\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 36\u001b[0;31m \u001b[0max\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mplot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mm\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     37\u001b[0m \u001b[0max\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mgrid\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     38\u001b[0m \u001b[0max\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mset_xlim\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m500\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mAttributeError\u001b[0m: 'numpy.ndarray' object has no attribute 'plot'"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA7AAAAD8CAYAAABOz0hIAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAby0lEQVR4nO3dX4ild30/8Pfnt+tCTUWtWa1sYk3L1jQFA3Ga2GrVtGg3obIIXiQVBRGWtKaUXhRDL+xFb1q8KbbRsEgQLzQX1egW1CiUVmmadmdLzD8b2UZrhhWyUVGqpWH187s4J3Ayzuw8s2dmzjw5rxccdp7n+T5nP/Nhhjefc545T3V3AAAAYL/7f4suAAAAAIYwwAIAADAKBlgAAABGwQALAADAKBhgAQAAGAUDLAAAAKOw5QBbVXdX1VNV9cgmx6uqPlxVZ6vqoaq6bubYsap6fHrsjp0sHACWlWwGYFkNeQf240mOXeT4TUmOTh8nknw0SarqQJI7p8evSXJrVV0zT7EAQBLZDMCS2nKA7e6vJPneRZYcT/KJnnggyUuq6pVJrk9ytruf6O5nktwzXQsAzEE2A7CsDu7AcxxJ8uTM9tp030b7b9jsSarqRCavEueyyy573dVXX70DpQGw7M6cOfN0dx9edB17TDYDsG/Nk807McDWBvv6Ivs31N0nk5xMkpWVlV5dXd2B0gBYdlX134uuYQFkMwD71jzZvBMD7FqSK2e2r0hyLsmhTfYDALtLNgPwvLQTt9E5leQ90088fH2SH3T3d5KcTnK0qq6qqkNJbpmuBQB2l2wG4Hlpy3dgq+pTSd6S5PKqWkvyF0lekCTdfVeSzye5OcnZJD9O8t7psQtVdXuS+5IcSHJ3dz+6C98DACwV2QzAstpygO3uW7c43knev8mxz2cSogDADpHNACyrnbiEGAAAAHadARYAAIBRMMACAAAwCgZYAAAARsEACwAAwCgYYAEAABgFAywAAACjYIAFAABgFAywAAAAjIIBFgAAgFEwwAIAADAKBlgAAABGwQALAADAKBhgAQAAGAUDLAAAAKNggAUAAGAUBg2wVXWsqh6vqrNVdccGx/+sqh6cPh6pqp9U1S9Mj32rqh6eHlvd6W8AAJaRbAZgGR3cakFVHUhyZ5K3JllLcrqqTnX3Y8+u6e4PJfnQdP3bk/xpd39v5mlu7O6nd7RyAFhSshmAZTXkHdjrk5zt7ie6+5kk9yQ5fpH1tyb51E4UBwBsSDYDsJSGDLBHkjw5s7023fczquqFSY4l+fTM7k7ypao6U1UnNvtPqupEVa1W1er58+cHlAUAS0s2A7CUhgywtcG+3mTt25P8y7pLlN7Q3dcluSnJ+6vqTRud2N0nu3ulu1cOHz48oCwAWFqyGYClNGSAXUty5cz2FUnObbL2lqy7RKm7z03/fSrJvZlc9gQAXDrZDMBSGjLAnk5ytKquqqpDmQThqfWLqurFSd6c5HMz+y6rqhc9+3WStyV5ZCcKB4AlJpsBWEpbfgpxd1+oqtuT3JfkQJK7u/vRqrptevyu6dJ3JPlSd/9o5vRXJLm3qp79vz7Z3V/cyW8AAJaNbAZgWVX3Zn8yszgrKyu9uuq2dADMr6rOdPfKousYO9kMwE6ZJ5uHXEIMAAAAC2eABQAAYBQMsAAAAIyCARYAAIBRMMACAAAwCgZYAAAARsEACwAAwCgYYAEAABgFAywAAACjYIAFAABgFAywAAAAjIIBFgAAgFEwwAIAADAKBlgAAABGwQALAADAKBhgAQAAGIVBA2xVHauqx6vqbFXdscHxt1TVD6rqwenjg0PPBQC2TzYDsIwObrWgqg4kuTPJW5OsJTldVae6+7F1S7/a3b9/iecCAAPJZgCW1ZB3YK9Pcra7n+juZ5Lck+T4wOef51wAYGOyGYClNGSAPZLkyZnttem+9X6zqr5WVV+oql/f5rmpqhNVtVpVq+fPnx9QFgAsLdkMwFIaMsDWBvt63fZ/JPml7r42yd8m+ew2zp3s7D7Z3SvdvXL48OEBZQHA0pLNACylIQPsWpIrZ7avSHJudkF3/7C7/2f69eeTvKCqLh9yLgCwbbIZgKU0ZIA9neRoVV1VVYeS3JLk1OyCqvrFqqrp19dPn/e7Q84FALZNNgOwlLb8FOLuvlBVtye5L8mBJHd396NVddv0+F1J3pnkD6vqQpL/TXJLd3eSDc/dpe8FAJaCbAZgWdUky/aXlZWVXl1dXXQZADwPVNWZ7l5ZdB1jJ5sB2CnzZPOQS4gBAABg4QywAAAAjIIBFgAAgFEwwAIAADAKBlgAAABGwQALAADAKBhgAQAAGAUDLAAAAKNggAUAAGAUDLAAAACMggEWAACAUTDAAgAAMAoGWAAAAEbBAAsAAMAoGGABAAAYhUEDbFUdq6rHq+psVd2xwfF3VdVD08f9VXXtzLFvVdXDVfVgVa3uZPEAsKxkMwDL6OBWC6rqQJI7k7w1yVqS01V1qrsfm1n2zSRv7u7vV9VNSU4muWHm+I3d/fQO1g0AS0s2A7CshrwDe32Ss939RHc/k+SeJMdnF3T3/d39/enmA0mu2NkyAYAZshmApTRkgD2S5MmZ7bXpvs28L8kXZrY7yZeq6kxVndjspKo6UVWrVbV6/vz5AWUBwNKSzQAspS0vIU5SG+zrDRdW3ZhJSL5xZvcbuvtcVb08yZer6j+7+ys/84TdJzO5vCkrKysbPj8AkEQ2A7CkhrwDu5bkypntK5KcW7+oql6b5GNJjnf3d5/d393npv8+leTeTC57AgAunWwGYCkNGWBPJzlaVVdV1aEktyQ5Nbugql6V5DNJ3t3d35jZf1lVvejZr5O8LckjO1U8ACwp2QzAUtryEuLuvlBVtye5L8mBJHd396NVddv0+F1JPpjkZUk+UlVJcqG7V5K8Ism9030Hk3yyu7+4K98JACwJ2QzAsqru/fcnLSsrK7266rZ0AMyvqs5MBzfmIJsB2CnzZPOQS4gBAABg4QywAAAAjIIBFgAAgFEwwAIAADAKBlgAAABGwQALAADAKBhgAQAAGAUDLAAAAKNggAUAAGAUDLAAAACMggEWAACAUTDAAgAAMAoGWAAAAEbBAAsAAMAoGGABAAAYBQMsAAAAozBogK2qY1X1eFWdrao7NjheVfXh6fGHquq6oecCANsnmwFYRlsOsFV1IMmdSW5Kck2SW6vqmnXLbkpydPo4keSj2zgXANgG2QzAshryDuz1Sc529xPd/UySe5IcX7fmeJJP9MQDSV5SVa8ceC4AsD2yGYCldHDAmiNJnpzZXktyw4A1RwaemySpqhOZvEKcJP9XVY8MqI3NXZ7k6UUXMXJ6OD89nJ8ezu81iy5gF8jmcfL7PD89nJ8ezk8P53fJ2TxkgK0N9vXANUPOnezsPpnkZJJU1Wp3rwyojU3o4fz0cH56OD89nF9VrS66hl0gm0dID+enh/PTw/np4fzmyeYhA+xakitntq9Icm7gmkMDzgUAtkc2A7CUhvwN7OkkR6vqqqo6lOSWJKfWrTmV5D3TTzx8fZIfdPd3Bp4LAGyPbAZgKW35Dmx3X6iq25Pcl+RAkru7+9Gqum16/K4kn09yc5KzSX6c5L0XO3dAXScv5ZvhOfRwfno4Pz2cnx7O73nXQ9k8Wno4Pz2cnx7OTw/nd8k9rO4N/+wFAAAA9pUhlxADAADAwhlgAQAAGIWFDbBVdayqHq+qs1V1xwbHq6o+PD3+UFVdt4g697MBPXzXtHcPVdX9VXXtIurc77bq48y636iqn1TVO/eyvjEY0sOqektVPVhVj1bVP+91jfvdgN/nF1fVP1TV16Y9fO8i6tzPquruqnpqs3uVypWtyeb5yeadIZvnJ5vnJ5vns2u53N17/sjkQyP+K8kvZ/Jx/l9Lcs26NTcn+UIm96t7fZJ/W0St+/UxsIe/leSl069v0sNL6+PMun/M5ENR3rnouvfTY+DP4kuSPJbkVdPtly+67v30GNjDP0/y19OvDyf5XpJDi659Pz2SvCnJdUke2eS4XLl4/2Tz3vRQNu9AH2fWyeZL7KFs3pEeyuaL93BXcnlR78Ben+Rsdz/R3c8kuSfJ8XVrjif5RE88kOQlVfXKvS50H9uyh919f3d/f7r5QCb3+uO5hvwsJskfJ/l0kqf2sriRGNLDP0jyme7+dpJ0tz4+15AedpIXVVUl+flMQvLC3pa5v3X3VzLpy2bkysXJ5vnJ5p0hm+cnm+cnm+e0W7m8qAH2SJInZ7bXpvu2u2aZbbc/78vkFQ6ea8s+VtWRJO9Ictce1jUmQ34WfzXJS6vqn6rqTFW9Z8+qG4chPfy7JL+W5FySh5P8SXf/dG/Ke96QKxcnm+cnm3eGbJ6fbJ6fbN59l5QpW94HdpfUBvvW389nyJplNrg/VXVjJiH5xl2taJyG9PFvknygu38yeYGNdYb08GCS1yX53SQ/l+Rfq+qB7v7Gbhc3EkN6+HtJHkzyO0l+JcmXq+qr3f3DXa7t+USuXJxsnp9s3hmyeX6yeX6yefddUqYsaoBdS3LlzPYVmbxysd01y2xQf6rqtUk+luSm7v7uHtU2JkP6uJLknmlAXp7k5qq60N2f3ZMK97+hv89Pd/ePkvyoqr6S5NokQnJiSA/fm+SvevJHI2er6ptJrk7y73tT4vOCXLk42Tw/2bwzZPP8ZPP8ZPPuu6RMWdQlxKeTHK2qq6rqUJJbkpxat+ZUkvdMP53q9Ul+0N3f2etC97Ete1hVr0rymSTv9mraprbsY3df1d2v7u5XJ/n7JH8kIJ9jyO/z55L8dlUdrKoXJrkhydf3uM79bEgPv53Jq+SpqlckeU2SJ/a0yvGTKxcnm+cnm3eGbJ6fbJ6fbN59l5QpC3kHtrsvVNXtSe7L5BO+7u7uR6vqtunxuzL5RLmbk5xN8uNMXuFgamAPP5jkZUk+Mn2F8kJ3ryyq5v1oYB+5iCE97O6vV9UXkzyU5KdJPtbdG36k+jIa+HP4l0k+XlUPZ3LJzQe6++mFFb0PVdWnkrwlyeVVtZbkL5K8IJErQ8jm+cnmnSGb5yeb5yeb57dbuVyTd7wBAABgf9vyEuKa4wa0NfAm1ADAcLIZgGU15G9gP57k2EWO35Tk6PRxIslHk6SqDiS5c3r8miS3VtU18xQLACSRzQAsqS0H2DluQDv0JtQAwDbIZgCW1U58iNNmN6DdaP8Nmz1JVZ3I5FXiXHbZZa+7+uqrd6A0AJbdmTNnnu7uw4uuY4/JZgD2rXmyeScG2M1uQLutG9N298kkJ5NkZWWlV1dXd6A0AJZdVf33omtYANkMwL41TzbvxAC72Q1oD22yHwDYXbIZgOelIR/itJXNbkA75Oa/AMDOk80APC9t+Q7spd6AdrOb/+7C9wAAS0U2A7Csthxgu/vWLY53kvdvcuzzmYQoALBDZDMAy2onLiEGAACAXWeABQAAYBQMsAAAAIyCARYAAIBRMMACAAAwCgZYAAAARsEACwAAwCgYYAEAABgFAywAAACjYIAFAABgFAywAAAAjIIBFgAAgFEwwAIAADAKBlgAAABGwQALAADAKBhgAQAAGIVBA2xVHauqx6vqbFXdscHxP6uqB6ePR6rqJ1X1C9Nj36qqh6fHVnf6GwCAZSSbAVhGB7daUFUHktyZ5K1J1pKcrqpT3f3Ys2u6+0NJPjRd//Ykf9rd35t5mhu7++kdrRwAlpRsBmBZDXkH9vokZ7v7ie5+Jsk9SY5fZP2tST61E8UBABuSzQAspSED7JEkT85sr033/YyqemGSY0k+PbO7k3ypqs5U1YnN/pOqOlFVq1W1ev78+QFlAcDSks0ALKUhA2xtsK83Wfv2JP+y7hKlN3T3dUluSvL+qnrTRid298nuXunulcOHDw8oCwCWlmwGYCkNGWDXklw5s31FknObrL0l6y5R6u5z03+fSnJvJpc9AQCXTjYDsJSGDLCnkxytqquq6lAmQXhq/aKqenGSNyf53My+y6rqRc9+neRtSR7ZicIBYInJZgCW0pafQtzdF6rq9iT3JTmQ5O7ufrSqbpsev2u69B1JvtTdP5o5/RVJ7q2qZ/+vT3b3F3fyGwCAZSObAVhW1b3Zn8wszsrKSq+uui0dAPOrqjPdvbLoOsZONgOwU+bJ5iGXEAMAAMDCGWABAAAYBQMsAAAAo2CABQAAYBQMsAAAAIyCARYAAIBRMMACAAAwCgZYAAAARsEACwAAwCgYYAEAABgFAywAAACjYIAFAABgFAywAAAAjIIBFgAAgFEwwAIAADAKgwbYqjpWVY9X1dmqumOD42+pqh9U1YPTxweHngsAbJ9sBmAZHdxqQVUdSHJnkrcmWUtyuqpOdfdj65Z+tbt//xLPBQAGks0ALKsh78Ben+Rsdz/R3c8kuSfJ8YHPP8+5AMDGZDMAS2nIAHskyZMz22vTfev9ZlV9raq+UFW/vs1zU1Unqmq1qlbPnz8/oCwAWFqyGYClNGSArQ329brt/0jyS919bZK/TfLZbZw72dl9srtXunvl8OHDA8oCgKUlmwFYSkMG2LUkV85sX5Hk3OyC7v5hd//P9OvPJ3lBVV0+5FwAYNtkMwBLacgAezrJ0aq6qqoOJbklyanZBVX1i1VV06+vnz7vd4ecCwBsm2wGYClt+SnE3X2hqm5Pcl+SA0nu7u5Hq+q26fG7krwzyR9W1YUk/5vklu7uJBueu0vfCwAsBdkMwLKqSZbtLysrK726urroMgB4HqiqM929sug6xk42A7BT5snmIZcQAwAAwMIZYAEAABgFAywAAACjYIAFAABgFAywAAAAjIIBFgAAgFEwwAIAADAKBlgAAABGwQALAADAKBhgAQAAGAUDLAAAAKNggAUAAGAUDLAAAACMggEWAACAUTDAAgAAMAoGWAAAAEZh0ABbVceq6vGqOltVd2xw/F1V9dD0cX9VXTtz7FtV9XBVPVhVqztZPAAsK9kMwDI6uNWCqjqQ5M4kb02yluR0VZ3q7sdmln0zyZu7+/tVdVOSk0lumDl+Y3c/vYN1A8DSks0ALKsh78Ben+Rsdz/R3c8kuSfJ8dkF3X1/d39/uvlAkit2tkwAYIZsBmApDRlgjyR5cmZ7bbpvM+9L8oWZ7U7ypao6U1UnNjupqk5U1WpVrZ4/f35AWQCwtGQzAEtpy0uIk9QG+3rDhVU3ZhKSb5zZ/YbuPldVL0/y5ar6z+7+ys88YffJTC5vysrKyobPDwAkkc0ALKkh78CuJblyZvuKJOfWL6qq1yb5WJLj3f3dZ/d397npv08luTeTy54AgEsnmwFYSkMG2NNJjlbVVVV1KMktSU7NLqiqVyX5TJJ3d/c3ZvZfVlUvevbrJG9L8shOFQ8AS0o2A7CUtryEuLsvVNXtSe5LciDJ3d39aFXdNj1+V5IPJnlZko9UVZJc6O6VJK9Icu9038Ekn+zuL+7KdwIAS0I2A7Csqnv//UnLyspKr666LR0A86uqM9PBjTnIZgB2yjzZPOQSYgAAAFg4AywAAACjYIAFAABgFAywAAAAjIIBFgAAgFEwwAIAADAKBlgAAABGwQALAADAKBhgAQAAGAUDLAAAAKNggAUAAGAUDLAAAACMggEWAACAUTDAAgAAMAoGWAAAAEbBAAsAAMAoDBpgq+pYVT1eVWer6o4NjldVfXh6/KGqum7ouQDA9slmAJbRlgNsVR1IcmeSm5Jck+TWqrpm3bKbkhydPk4k+eg2zgUAtkE2A7CshrwDe32Ss939RHc/k+SeJMfXrTme5BM98UCSl1TVKweeCwBsj2wGYCkdHLDmSJInZ7bXktwwYM2RgecmSarqRCavECfJ/1XVIwNqY3OXJ3l60UWMnB7OTw/np4fze82iC9gFsnmc/D7PTw/np4fz08P5XXI2Dxlga4N9PXDNkHMnO7tPJjmZJFW12t0rA2pjE3o4Pz2cnx7OTw/nV1Wri65hF8jmEdLD+enh/PRwfno4v3myecgAu5bkypntK5KcG7jm0IBzAYDtkc0ALKUhfwN7OsnRqrqqqg4luSXJqXVrTiV5z/QTD1+f5Afd/Z2B5wIA2yObAVhKW74D290Xqur2JPclOZDk7u5+tKpumx6/K8nnk9yc5GySHyd578XOHVDXyUv5ZngOPZyfHs5PD+enh/N73vVQNo+WHs5PD+enh/PTw/ldcg+re8M/ewEAAIB9ZcglxAAAALBwBlgAAABGYWEDbFUdq6rHq+psVd2xwfGqqg9Pjz9UVdctos79bEAP3zXt3UNVdX9VXbuIOve7rfo4s+43quonVfXOvaxvDIb0sKreUlUPVtWjVfXPe13jfjfg9/nFVfUPVfW1aQ/fu4g697OquruqntrsXqVyZWuyeX6yeWfI5vnJ5vnJ5vnsWi53954/MvnQiP9K8suZfJz/15Jcs27NzUm+kMn96l6f5N8WUet+fQzs4W8leen065v08NL6OLPuHzP5UJR3Lrru/fQY+LP4kiSPJXnVdPvli657Pz0G9vDPk/z19OvDSb6X5NCia99PjyRvSnJdkkc2OS5XLt4/2bw3PZTNO9DHmXWy+RJ7KJt3pIey+eI93JVcXtQ7sNcnOdvdT3T3M0nuSXJ83ZrjST7REw8keUlVvXKvC93Htuxhd9/f3d+fbj6Qyb3+eK4hP4tJ8sdJPp3kqb0sbiSG9PAPknymu7+dJN2tj881pIed5EVVVUl+PpOQvLC3Ze5v3f2VTPqyGblycbJ5frJ5Z8jm+cnm+cnmOe1WLi9qgD2S5MmZ7bXpvu2uWWbb7c/7MnmFg+faso9VdSTJO5LctYd1jcmQn8VfTfLSqvqnqjpTVe/Zs+rGYUgP/y7JryU5l+ThJH/S3T/dm/KeN+TKxcnm+cnmnSGb5yeb5yebd98lZcqW94HdJbXBvvX38xmyZpkN7k9V3ZhJSL5xVysapyF9/JskH+jun0xeYGOdIT08mOR1SX43yc8l+deqeqC7v7HbxY3EkB7+XpIHk/xOkl9J8uWq+mp3/3CXa3s+kSsXJ5vnJ5t3hmyen2yen2zefZeUKYsaYNeSXDmzfUUmr1xsd80yG9Sfqnptko8luam7v7tHtY3JkD6uJLlnGpCXJ7m5qi5092f3pML9b+jv89Pd/aMkP6qqryS5NomQnBjSw/cm+aue/NHI2ar6ZpKrk/z73pT4vCBXLk42z0827wzZPD/ZPD/ZvPsuKVMWdQnx6SRHq+qqqjqU5JYkp9atOZXkPdNPp3p9kh9093f2utB9bMseVtWrknwmybu9mrapLfvY3Vd196u7+9VJ/j7JHwnI5xjy+/y5JL9dVQer6oVJbkjy9T2ucz8b0sNvZ/IqearqFUlek+SJPa1y/OTKxcnm+cnmnSGb5yeb5yebd98lZcpC3oHt7gtVdXuS+zL5hK+7u/vRqrptevyuTD5R7uYkZ5P8OJNXOJga2MMPJnlZko9MX6G80N0ri6p5PxrYRy5iSA+7++tV9cUkDyX5aZKPdfeGH6m+jAb+HP5lko9X1cOZXHLzge5+emFF70NV9akkb0lyeVWtJfmLJC9I5MoQsnl+snlnyOb5yeb5yeb57VYu1+QdbwAAANjfFnUJMQAAAGyLARYAAIBRMMACAAAwCgZYAAAARsEACwAAwCgYYAEAABgFAywAAACj8P8Brzjsh4wuiCoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1152x288 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import math\n",
    "\n",
    "def getSin(amp,freq,phase,sampleList):\n",
    "    return amp*np.sin(-2*math.pi*freq*sampleList+phase)\n",
    "def getCos(amp,freq,phase,sampleList):\n",
    "    return amp*np.cos(-2*math.pi*freq*sampleList+phase)\n",
    "# 1. 获得混合波形\n",
    "srate=3000\n",
    "t=np.linspace(0,1,srate)\n",
    "s1=getSin(amp=5, freq=30, phase=0,sampleList=t)    \n",
    "s2=getCos(amp=3, freq=5, phase=0,sampleList=t)  \n",
    "s3=getSin(amp=10, freq=100, phase=0,sampleList=t)\n",
    "s4=getCos(amp=20, freq=120, phase=0,sampleList=t)  \n",
    "m=s1+s2+s3+s4\n",
    "# 2. 获得傅里叶系数\n",
    "fCoefs = np.fft.fft(m,srate)\n",
    "# 3. 获得振幅列表：每一个绕线的重心到原点的距离\n",
    "amp_list=2*np.abs(fCoefs/srate)\n",
    "# 把频率轴从0~1000 转变成 0~499 然后 -500~-1\n",
    "freqs = np.fft.fftfreq(len(amp_list), 1/srate)\n",
    "# 然后把 频率轴 和 数据 都变成 0hz 在中间，向左是负频率，向右是正频率的形式\n",
    "amp_shifted=np.fft.fftshift(amp_list)\n",
    "freq_shift=np.fft.fftshift(freqs)\n",
    "amp_shifted_copy = np.copy(amp_shifted)#获取图像\n",
    "amp_shifted_copy[amp_shifted < 0.5] = 0#使用掩模\n",
    "amp_shifted_copy[(freq_shift >= 110) | (freq_shift <= -110)] = 0\n",
    "#逆傅里叶变化\n",
    "amp_shifted_ishift=np.fft.ifftshift(amp_shifted_copy)\n",
    "amp_shifted3 = amp_shifted_ishift * srate / 2\n",
    "amp_shifted_ifft=np.fft.ifft(amp_shifted3)\n",
    "# 把振幅列表画出来\n",
    "fig,ax=plt.subplots(2,2,figsize=(16,4))\n",
    "#左上角子图\n",
    "ax[0].plot(m)\n",
    "ax[0].grid()\n",
    "ax[0].set_xlim([0,500])\n",
    "#右上角子图\n",
    "ax[1].plot(np.real(amp_shifted_ifft))\n",
    "ax[1].grid()\n",
    "ax[1].set_xlim([0,500])\n",
    "#左下角子图\n",
    "ax[2].grid()\n",
    "ax[2].stem(freq_shift,amp_shifted)\n",
    "ax[2].set_xlim([-150,150])\n",
    "#右下角子图\n",
    "ax[3].stem(freq_shift,amp_shifted_copy)\n",
    "ax[3].set_xlim([-150,150])\n",
    "ax[3].set_ylim([0, 20])\n",
    "ax[3].grid()\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7045e406-7ba7-4435-be05-3f22b791c24c",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "27053a99-22a9-43ba-a6fb-d1392f1d05e7",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
